#read global mean values from file
library(openxlsx)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
library(tidyr)

setwd("~/Documents/Promotion/RS-fPET/Analysen/flexible_factorial_multiframes")

#import globals
globals<-c(3319.94816708075,
           3418.38236634798,
           3497.18808493240,
           3577.02786408874,
           3659.49101254483,
           3756.67744453291,
           3809.87836042326,
           3898.57802230989,
           3999.31684429413,
           4067.43533060166,
           4141.67033527087,
           4227.20607622073,
           4302.37920166778,
           4387.44581454623,
           4470.22750970148,
           4523.49026427971,
           4592.58739233506,
           4638.80380850590,
           4657.24528367891,
           4721.00248055408,
           4752.43926965674,
           4810.47495914204,
           4811.32067331599,
           4859.34592704655,
           4898.67618278456,
           4943.01107473273,
           4949.49331189629,
           5007.50625781463,
           5060.08885486849,
           5074.08189170276,
           3463.21449367145,
           3555.81843063848,
           3639.67267486807,
           3734.10450296421,
           3801.91657696198,
           3895.99314360205,
           3988.07394164183,
           4041.65026866628,
           4110.65100663420,
           4192.70609698503,
           4292.16732758294,
           4377.10338527132,
           4501.34535697103,
           4567.50772308234,
           4678.20219486598,
           4740.08440999034,
           4811.25442556765,
           4950.06230568339,
           5021.62014786680,
           5072.87333014877,
           5129.44649611283,
           5167.32522983417,
           5239.80294998409,
           5266.96495412684,
           5340.24395329966,
           5332.45039552743,
           5368.89692345910,
           5442.17776526553,
           5495.03385877380,
           5531.23103340016,
           2724.70573079751,
           2801.11826425685,
           2841.94192562966,
           2895.85016620547,
           2979.67118069875,
           3037.26791561374,
           3084.16229164628,
           3141.66176872056,
           3213.61968135634,
           3251.18932015025,
           3354.87844066079,
           3415.80844849723,
           3442.50561038114,
           3522.19001883295,
           3616.54687438812,
           3650.15476273535,
           3676.84572754950,
           3726.55949529830,
           3758.13322674733,
           3766.31175750036,
           3811.14747796121,
           3847.24413594188,
           3875.38767129541,
           3916.72257360192,
           3935.29158905346,
           4008.14502065617,
           4007.56326332832,
           4017.61149758045,
           4064.51578663252,
           4080.07384297825,
           2757.73491821340,
           2843.03943303919,
           2901.70930179077,
           2947.41837937330,
           3023.03942607961,
           3082.33762490042,
           3153.83327096665,
           3196.95147980858,
           3247.36186139382,
           3296.83185193454,
           3335.99736248624,
           3439.28392699495,
           3502.18362923834,
           3559.72237660085,
           3601.11486783002,
           3671.49462353559,
           3753.15602501133,
           3816.38627929906,
           3895.21495904310,
           3943.57474063943,
           3989.94693071748,
           3997.53670790789,
           4018.32091037161,
           4053.80621386861,
           4087.91674402029,
           4123.19040303367,
           4147.52051605164,
           4187.17629939089,
           4203.47891858052,
           4236.66892088889,
           2764.84129842810,
           2831.63049094681,
           2868.23098580737,
           2943.11071953718,
           3006.72909143078,
           3050.72280082254,
           3128.98465889064,
           3191.23539233331,
           3225.09571096278,
           3298.38408161486,
           3361.11676816427,
           3431.46253842500,
           3487.51032283339,
           3553.70501444626,
           3626.84549881193,
           3684.91259194636,
           3729.49641294629,
           3819.25979116901,
           3857.36478212511,
           3901.48035701993,
           3944.02609392190,
           3979.57068712671,
           4011.23892713366,
           4037.22939138019,
           4072.33272491402,
           4098.10871243931,
           4119.79011274126,
           4136.17873738423,
           4179.31292774336,
           4202.42435953202,
           2723.64144111185,
           2773.95236388179,
           2833.74471530761,
           2898.16178133546,
           2949.09968331363,
           2992.62720962392,
           3063.71490928920,
           3109.14662896132,
           3196.59198556991,
           3250.88092906892,
           3313.01084222412,
           3375.63362340496,
           3447.15665147918,
           3499.09516315241,
           3569.49313945483,
           3629.23166233416,
           3692.81618154316,
           3751.45473922890,
           3838.66954179513,
           3866.62543219542,
           3933.03838373645,
           4022.04457899870,
           4092.28405984656,
           4147.63029732313,
           4216.44309019632,
           4260.31808141891,
           4302.52608464204,
           4308.82993789807,
           4370.55669470633,
           4422.24245858151,
           4579.32950466242,
           4695.54596611744,
           4826.28313527527,
           4941.78658796845,
           5070.41525124966,
           5204.68016355390,
           5312.88557906887,
           5414.32055575324,
           5552.69664562181,
           5665.88084524768,
           5770.85689811420,
           5910.12266769570,
           6055.05070816002,
           6183.28227574275,
           6321.77272297455,
           6453.90573384614,
           6567.16347121194,
           6693.51714455732,
           6845.42074371483,
           6978.83467842561,
           7147.19165375549,
           7267.04924369983,
           7397.61052857497,
           7550.31721321434,
           7682.23140108729,
           7817.06622310416,
           7914.70280347762,
           7989.47077624116,
           8079.99002331678,
           8132.24249384481,
           3056.76942796260,
           3079.33489472465,
           3069.81469325774,
           3102.69448316216,
           3114.07162722644,
           3123.54793197637,
           3130.63511890437,
           3131.36310117853,
           3162.91902240081,
           3145.20460355814,
           3165.47996677631,
           3180.30155220493,
           3193.02739021486,
           3193.88239047159,
           3214.11538545042,
           3210.46038037488,
           3232.77169162684,
           3235.82289734160,
           3252.38705476565,
           3260.45284295437,
           3259.55073254645,
           3272.28148583680,
           3290.35088447513,
           3301.16464624867,
           3309.35639926200,
           3309.53466245413,
           3318.78537516829,
           3334.99442146378,
           3334.28818751421,
           3350.24857297137,
           2700.62312867496,
           2755.06339996951,
           2834.28807717885,
           2884.48202811316,
           2979.37882914823,
           3031.58182553708,
           3102.96287881290,
           3159.86628633659,
           3212.49067008733,
           3313.03134998421,
           3393.99219187691,
           3456.27230331572,
           3533.04683923934,
           3603.42854410941,
           3673.28194057120,
           3767.63311059231,
           3850.06582376081,
           3913.29520744008,
           4010.54489933994,
           4080.26492847487,
           4159.40926760452,
           4242.93840701548,
           4304.85606515297,
           4396.46639206363,
           4490.70490701455,
           4563.65927237689,
           4626.59600565548,
           4727.00286163860,
           4793.77030095492,
           4881.22658436441,
           4672.95514577808,
           4783.40033970401,
           4886.54963545112,
           5000.74606887873,
           5136.50447848172,
           5256.37723139465,
           5356.75318282963,
           5486.23087647404,
           5597.94950718557,
           5734.09424436729,
           5831.52298771340,
           5945.27164279603,
           6070.81041059454,
           6180.51836916879,
           6320.35049299492,
           6444.03620962951,
           6532.63823219008,
           6674.69319906543,
           6806.41546297530,
           6870.46406698631,
           7032.25620628118,
           7172.34339203785,
           7273.25446131122,
           7401.24180403076,
           7502.58459165665,
           7684.58555921490,
           7775.45281498599,
           7920.94821575943,
           8049.51885025144,
           8185.58483663639,
           2095.06239833971,
           2170.62447948463,
           2197.44574060888,
           2224.80754885885,
           2284.58267989466,
           2340.99848651292,
           2390.47434655120,
           2439.74748694997,
           2454.07915972278,
           2492.87999071611,
           2484.68442805821,
           2582.48439599537,
           2644.65058151161,
           2713.51234752882,
           2752.70857391481,
           2835.28138078382,
           2849.18176651456,
           2820.99475163803,
           2885.57492724856,
           2980.16734319188,
           3040.64176818402,
           3096.62513823249,
           3131.98450663043,
           3186.28623348476,
           3227.99272130983,
           3282.51509868892,
           3398.71277817170,
           3425.58471209891,
           3432.67547188267,
           3551.87896379541,
           3257.50458367897,
           3317.40485680203,
           3365.07081818868,
           3447.94980640515,
           3549.27168614472,
           3596.00705635560,
           3717.30031242694,
           3802.13508250150,
           3917.76123578062,
           3960.42220524238,
           4063.21324091555,
           4131.31270674617,
           4162.88910898933,
           4202.50848785036,
           4252.12673433285,
           4287.91308468453,
           4312.40377621723,
           4368.89807513131,
           4384.84724778483,
           4432.64163500377,
           4504.28195761995,
           4528.48860036807,
           4587.52904059029,
           4665.55910797011,
           4770.64986654619,
           4846.35881427521,
           4930.32911283703,
           5039.10203474194,
           5049.84798240388,
           5087.07554471904,
           3100.84925019087,
           3179.93087804215,
           3254.30661126453,
           3336.54711194584,
           3412.62460927669,
           3495.78105421507,
           3579.98975234024,
           3654.23387396665,
           3726.73680647839,
           3817.79956559307,
           3900.39001548916,
           3978.66007161610,
           4058.09336105587,
           4155.00022625346,
           4229.02555258606,
           4327.28173478395,
           4398.08719229951,
           4497.73414778229,
           4572.18681449414,
           4660.38519731358,
           4769.72161316680,
           4852.92833046445,
           4945.37497241966,
           5019.09043501711,
           5105.12595091168,
           5191.87801954082,
           5294.72453935502,
           5418.08700713833,
           5438.25379113736,
           5506.43308116922,
           2479.84633469280,
           2558.45863745364,
           2586.74015175475,
           2640.22563135263,
           2677.61654918587,
           2684.53495994127,
           2761.02425179149,
           2829.77632823087,
           2894.06121548796,
           2959.52543335916,
           3002.91526262532,
           3065.29641932065,
           3140.64704133912,
           3185.89986248209,
           3255.77755388229,
           3297.32204327621,
           3398.19865565633,
           3436.13304742216,
           3498.99201231784,
           3579.74691009578,
           3675.43046536305,
           3747.11089701406,
           3818.35303221822,
           3894.58254220171,
           3872.62796097249,
           3916.38682537284,
           3958.77267715184,
           4041.37263780339,
           4102.37920044991,
           4136.85644470996,
           4441.95439767401,
           4534.30210515798,
           4678.90025586019,
           4767.49654857554,
           4890.32460045874,
           5003.42653992291,
           5134.03690854965,
           5256.05442907187,
           5376.58332146329,
           5497.88288118729,
           5618.94434752291,
           5732.84776340330,
           5848.96387339833,
           5980.89973344712,
           6083.74199874933,
           6211.08955242319,
           6346.15326389427,
           6474.02361414170,
           6578.85093397235,
           6703.96094708056,
           6861.13603246712,
           6953.17563287754,
           7093.53304435392,
           7205.89217918149,
           7344.35016682937,
           7471.63393466732,
           7590.95893123734,
           7729.12903252164,
           7857.34999368733,
           7943.76762260336,
           2845.54900754000,
           2917.05861909940,
           2987.82080796886,
           3088.78939026572,
           3166.29768659137,
           3294.31516871021,
           3364.41859549549,
           3424.46823041624,
           3518.04184246446,
           3589.37473511420,
           3671.54886693156,
           3708.77400784638,
           3854.34305464441,
           3955.07087503823,
           4015.42241716590,
           4093.24139781126,
           4186.00609458063,
           4247.02843295876,
           4356.07447440638,
           4474.60110418588,
           4524.00969304423,
           4610.55937407860,
           4690.56816891014,
           4769.46563289049,
           4883.92801080927,
           4941.39255318814,
           5045.64882476297,
           5139.80823816589,
           5209.26695196210,
           5279.03926299111,
           4292.51739610408,
           4389.49590255018,
           4501.81899080335,
           4618.92418680273,
           4723.12900484050,
           4814.08011313969,
           4928.43685016322,
           5009.06197088842,
           5109.35900078647,
           5218.30248947565,
           5321.69968506417,
           5433.34009990898,
           5555.47059522491,
           5659.89446570111,
           5753.96834632869,
           5872.75166724424,
           5963.51080422932,
           6067.71432420295,
           6174.02272654688,
           6282.47734960793,
           6391.03490259500,
           6478.86162117558,
           6629.30818143491,
           6722.97071202891,
           6826.83871857741,
           6939.83085226035,
           6982.91069179954,
           7068.08218962183,
           7151.68976043225,
           7193.55096766104,
           5887.87741559224,
           6056.86700499497,
           6197.39648564736,
           6347.75534685844,
           6469.03344320978,
           6626.83431506232,
           6746.21022048908,
           6911.43442595896,
           7068.51452351016,
           7213.45139808186,
           7386.82150453702,
           7515.36539326744,
           7662.24517212888,
           7849.03950296941,
           7979.44192855476,
           8178.69989581313,
           8344.61598216003,
           8472.14063858022,
           8647.75666381938,
           8777.26311492773,
           8961.04365313878,
           9099.22656682185,
           9280.20790468723,
           9435.01200524830,
           9615.09029284948,
           9754.41708808525,
           9918.30715733493,
           10087.8871693458,
           10233.5587514870,
           10396.0939874438,
           4579.33252094656,
           4695.59838833431,
           4826.31285687290,
           4941.81619588835,
           5070.39000863408,
           5204.71286680576,
           5312.91851661789,
           5414.29327049736,
           5552.70054926923,
           5665.85228663364,
           5770.89415883801,
           5910.12671945010,
           6054.95325438077,
           6183.35557584894,
           6321.77599960135,
           6453.94710444522,
           6567.24369094516,
           6693.59671744208,
           6845.39040265081,
           6978.79944364394,
           7147.15838049054,
           7267.01878559355,
           7397.61876803269,
           7550.28485210578,
           7682.28374193943,
           7817.02954160814,
           7914.70926271497,
           7989.47730817632,
           8079.95198319183,
           8132.29674717295,
           2616.55968027273,
           2698.68085193695,
           2744.80825447971,
           2807.21240724668,
           2869.56236155191,
           2935.36956747497,
           3003.76290041638,
           3060.36600456702,
           3146.01010207412,
           3201.57355448887,
           3253.39179866598,
           3340.17771251816,
           3401.46047397038,
           3478.61845134193,
           3546.64775305503,
           3607.64234635013,
           3642.99560786268,
           3729.26997034677,
           3808.84973295660,
           3841.80296899294,
           3942.10504904216,
           4000.50626008020,
           4054.55306936458,
           4128.74750839528,
           4196.84634100951,
           4245.88514064596,
           4305.38652995404,
           4385.65969013307,
           4417.00809879208,
           4478.64807433623,
           2910.08967780965,
           2990.63038056923,
           3042.83788963236,
           3128.31449581930,
           3199.56234745909,
           3279.30696810227,
           3365.80922260096,
           3410.87173234437,
           3513.94324672892,
           3574.90224974548,
           3636.05443030603,
           3721.34988897654,
           3789.69704889620,
           3853.17000109644,
           3946.32729302398,
           4008.27140365347,
           4093.37674560777,
           4171.94404835394,
           4234.02814473147,
           4322.59180003881,
           4421.90279467229,
           4486.63703164228,
           4562.90335489862,
           4672.59998682407,
           4731.36916801668,
           4822.48990138430,
           4883.83057986733,
           5004.51795504542,
           5071.95059849855,
           5165.68930788593,
           1595.02197722526,
           1638.82506120646,
           1662.39544284796,
           1704.78425724697,
           1740.44876228869,
           1780.07226154444,
           1813.72924692219,
           1856.26259312770,
           1906.68252340686,
           1932.43421617007,
           1967.22414590665,
           2030.81182339258,
           2053.92154281495,
           2101.58751286084,
           2136.26950565587,
           2185.34959929394,
           2239.77938375152,
           2238.16026894445,
           2290.93066517181,
           2336.21477256619,
           2370.54632018562,
           2431.34782680987,
           2461.42970551090,
           2463.31816783571,
           2526.64268090659,
           2576.27909282674,
           2606.27789975522,
           2657.63074971418,
           2710.45023055670,
           2728.84888161173,
           2830.38874262949,
           2922.19570724204,
           2972.47775759304,
           3052.22312356074,
           3091.90490437638,
           3169.83574380782,
           3230.84049794811,
           3281.36689162791,
           3352.91744841731,
           3407.36221083120,
           3481.71626080718,
           3534.25284706631,
           3589.16389475647,
           3664.05356446221,
           3714.72174200627,
           3770.83051925415,
           3859.09081227308,
           3950.61670148412,
           4025.95039420709,
           4099.30461256301,
           4198.73201656123,
           4264.35058088041,
           4334.36027917912,
           4416.65824384271,
           4456.74174673857,
           4532.41071578495,
           4636.55247795948,
           4683.28337572792,
           4751.12525029651,
           4786.72295484841,
           2616.90226029821,
           2687.99153768661,
           2765.20163477295,
           2825.36907028436,
           2888.27903816430,
           2952.15180143661,
           3027.93972279379,
           3097.90551229701,
           3180.62563989970,
           3241.68459941748,
           3310.43813694709,
           3387.71914287298,
           3429.33094047800,
           3527.34381768486,
           3600.61469575232,
           3658.24601944496,
           3740.16446202610,
           3805.26611363195,
           3867.99829405987,
           3919.16869566141,
           4012.12919413232,
           4075.38318184229,
           4153.24682474383,
           4209.56464903985,
           4281.06435434925,
           4368.61650095518,
           4441.33064997784,
           4511.04393370240,
           4591.35236671113,
           4660.41253262556,
           3401.74954024533,
           3470.92629357341,
           3546.04232399719,
           3658.75923587124,
           3728.84040983163,
           3808.00988761300,
           3908.16113098115,
           3988.31563402157,
           4086.06295940815,
           4178.56025872323,
           4263.53224166557,
           4342.00505681526,
           4433.13288197271,
           4512.91342570933,
           4600.93605821228,
           4671.44885840151,
           4797.44882853567,
           4873.32939619339,
           4984.70755605811,
           5051.13892464322,
           5122.22761752905,
           5222.70992542209,
           5300.63075602789,
           5428.82151907296,
           5528.33327599924,
           5595.13578139300,
           5689.25021625255,
           5766.28779112918,
           5861.11393302149,
           5926.53774068577,
           3222.59646860021,
           3338.35298544917,
           3408.64082530093,
           3475.53994809764,
           3596.05833640940,
           3671.70707456757,
           3744.54674953722,
           3857.89743163596,
           3961.39389545491,
           4050.15292937940,
           4132.45286845147,
           4201.15240970621,
           4276.25891870123,
           4352.26611336383,
           4477.30659872174,
           4564.13647048005,
           4648.82179053384,
           4707.84952303158,
           4827.00996343507,
           4888.82586608797,
           4989.13281016028,
           5052.23910247192,
           5172.91218690437,
           5269.21922942399,
           5346.00625930406,
           5436.43195410137,
           5545.65061055734,
           5595.30060931343,
           5700.27646432542,
           5765.46523956464,
           2457.66540424913,
           2526.78924200534,
           2597.57297735770,
           2651.42326132279,
           2709.52522562361,
           2764.50236610288,
           2828.78536941098,
           2887.30069481355,
           2955.35541733852,
           3042.40184578059,
           3108.94890230372,
           3170.88721522308,
           3239.38520674147,
           3312.77400626023,
           3385.02818353278,
           3486.98273682513,
           3567.72053557674,
           3655.75021520328,
           3748.92138712446,
           3753.46501947168,
           3854.78289984392,
           3875.25529380633,
           3971.53510476111,
           4031.87521561430,
           4103.70865248656,
           4170.53428709115,
           4233.28388026583,
           4309.40304256763,
           4398.00708602172,
           4406.17706750327)

#read data
FlexData<-read.xlsx(xlsxFile = "flex_fact_new2025_02_05.xlsx", sheet = "Tabelle1",
                           colNames = TRUE)

#devide values element-wise by global mean
gmFlexData<-FlexData/globals


#define subject variable subject id 
subject_id<-c(rep("RS_fPET001",30), rep("RS_fPET002",30), rep("RS_fPET003",30),rep("RS_fPET004",30), rep("RS_fPET005",30),rep("RS_fPET006",30),rep("RS_fPET007",30),rep("RS_fPET008",30),rep("RS_fPET010",30),
              rep("RS_fPET011",30),rep("RS_fPET012",30),rep("RS_fPET014",30),rep("RS_fPET015",30),rep("RS_fPET016",30),
              rep("RS_fPET017",30),
              rep("RS_fPET018",30),rep("RS_fPET019",30),rep("RS_fPET020",30), rep("RS_fPET021",30),rep("RS_fPET022",30),rep("RS_fPET023",30),rep("RS_fPET024",30),rep("RS_fPET025",30),
              rep("RS_fPET026",30),rep("RS_fPET027",30),rep("RS_fPET028",30),rep("RS_fPET029",30))


#add and convert record id into factor
gmFlexData$record_id<-subject_id
gmFlexData$record_id<-factor(gmFlexData$record_id)

#add diagnose to_id gm removed data set 
gmFlexData$diagnose[gmFlexData$record_id=="RS_fPET001"|gmFlexData$record_id=="RS_fPET002"|gmFlexData$record_id=="RS_fPET003"|gmFlexData$record_id=="RS_fPET005"|gmFlexData$record_id=="RS_fPET007"|gmFlexData$record_id=="RS_fPET008"|gmFlexData$record_id=="RS_fPET010"|gmFlexData$record_id=="RS_fPET014"|gmFlexData$record_id=="RS_fPET016"|gmFlexData$record_id=="RS_fPET022"|gmFlexData$record_id=="RS_fPET025"|gmFlexData$record_id=="RS_fPET026"|gmFlexData$record_id=="RS_fPET028"|gmFlexData$record_id=="RS_fPET029"]= "PD"

gmFlexData$diagnose[gmFlexData$record_id=="RS_fPET004"|gmFlexData$record_id=="RS_fPET006"|gmFlexData$record_id=="RS_fPET011"|gmFlexData$record_id=="RS_fPET012"|gmFlexData$record_id=="RS_fPET015"|gmFlexData$record_id=="RS_fPET017"|gmFlexData$record_id=="RS_fPET018"|gmFlexData$record_id=="RS_fPET019"|gmFlexData$record_id=="RS_fPET020"|gmFlexData$record_id=="RS_fPET021"|gmFlexData$record_id=="RS_fPET023"|gmFlexData$record_id=="RS_fPET024"|gmFlexData$record_id=="RS_fPET027"]= "HC"

#convert diagnose into factor
gmFlexData$diagnose<-factor(gmFlexData$diagnose)
summary(gmFlexData$diagnose)
##  HC  PD 
## 390 420
#add timepoints
time<-c(rep(61:90, 27)) #repeat time points 61- 90 27 x
gmFlexData$time<-time
#save activity-time plots for all subjects with mean line per group 
for (i in 1:ncol(gmFlexData[,1:12])) {
  MeanTime <- aggregate(gmFlexData[,i] ~ time + diagnose, data = gmFlexData, FUN = mean)
  MeanTime_df <- cbind(MeanTime[-ncol(MeanTime)], MeanTime[[ncol(MeanTime)]])
  colnames(MeanTime_df)<-c("time","diagnose",colnames(gmFlexData)[i])
  for (j in 3:ncol(MeanTime_df)){
    plot <- ggplot(gmFlexData, aes(x = time, y = gmFlexData[,i])) +
      geom_line(aes(color = diagnose, group = record_id), linetype = 3) +
     
      geom_line(data = MeanTime_df, aes(x = time, y = MeanTime_df[,j], color = diagnose), linewidth = 1) + 
     
        theme_bw()+
      labs(title = "", x = "Time (Min)",
           y = expression("Relative ["^18*"F]-FDG Uptake"))+
      scale_color_manual(values=c("#1F78B4",  "#BFD200"))+
        theme(legend.position = c(0.1,0.85))+
       theme(legend.title = element_blank())+
         theme(text = element_text(size = 25),
        axis.text.x = element_text(size=20),
        axis.text.y = element_text(size=20))+
    show(plot)
    ggsave(plot,file=paste0("gmFlexData",colnames(gmFlexData)[i],".png"))}}
## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## function (x, y, ...) 
## UseMethod("plot")
## <bytecode: 0x127c88c40>
## <environment: namespace:base>
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

setwd(“~/Documents/Promotion/RS-fPET”)

#change wd


#mean uptake per subject 
MeanUptake <- aggregate(gmFlexData[,c(1:12)], by= list(subject=gmFlexData$record_id), FUN = mean)

#add diagnose to_id gm removed data set 
MeanUptake$diagnose[MeanUptake$subject=="RS_fPET001"|MeanUptake$subject=="RS_fPET002"|MeanUptake$subject=="RS_fPET003"|MeanUptake$subject=="RS_fPET005"|MeanUptake$subject=="RS_fPET007"|MeanUptake$subject=="RS_fPET008"|MeanUptake$subject=="RS_fPET010"|MeanUptake$subject=="RS_fPET014"|MeanUptake$subject=="RS_fPET016"|MeanUptake$subject=="RS_fPET022"|MeanUptake$subject=="RS_fPET025"|MeanUptake$subject=="RS_fPET026"|MeanUptake$subject=="RS_fPET028"|MeanUptake$subject=="RS_fPET029"]= "PD"

MeanUptake$diagnose[MeanUptake$subject=="RS_fPET004"|MeanUptake$subject=="RS_fPET006"|MeanUptake$subject=="RS_fPET011"|MeanUptake$subject=="RS_fPET012"|MeanUptake$subject=="RS_fPET015"|MeanUptake$subject=="RS_fPET017"|MeanUptake$subject=="RS_fPET018"|MeanUptake$subject=="RS_fPET019"|MeanUptake$subject=="RS_fPET020"|MeanUptake$subject=="RS_fPET021"|MeanUptake$subject=="RS_fPET023"|MeanUptake$subject=="RS_fPET024"|MeanUptake$subject=="RS_fPET027"]= "HC"


for (i in 2:ncol(MeanUptake[,2:14])){
boxplot(MeanUptake[,i]~diagnose, data=MeanUptake, main = paste0("Mean",colnames(MeanUptake)[i]))
results<-t.test(MeanUptake[,i]~diagnose, data=MeanUptake)
print(results)}

## 
##  Welch Two Sample t-test
## 
## data:  MeanUptake[, i] by diagnose
## t = 3.4257, df = 22.716, p-value = 0.002339
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
##  0.04290399 0.17393796
## sample estimates:
## mean in group HC mean in group PD 
##         1.492089         1.383668

## 
##  Welch Two Sample t-test
## 
## data:  MeanUptake[, i] by diagnose
## t = 2.2811, df = 21.302, p-value = 0.03293
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
##  0.01016617 0.21800837
## sample estimates:
## mean in group HC mean in group PD 
##         2.061555         1.947468

## 
##  Welch Two Sample t-test
## 
## data:  MeanUptake[, i] by diagnose
## t = 2.8208, df = 22.517, p-value = 0.009816
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
##  0.03491988 0.22785672
## sample estimates:
## mean in group HC mean in group PD 
##         1.736867         1.605478

## 
##  Welch Two Sample t-test
## 
## data:  MeanUptake[, i] by diagnose
## t = 3.0073, df = 22.5, p-value = 0.00638
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
##  0.04186964 0.22715419
## sample estimates:
## mean in group HC mean in group PD 
##         2.071889         1.937378

## 
##  Welch Two Sample t-test
## 
## data:  MeanUptake[, i] by diagnose
## t = 3.2926, df = 24.783, p-value = 0.002981
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
##  0.05223899 0.22694852
## sample estimates:
## mean in group HC mean in group PD 
##         2.007700         1.868106

## 
##  Welch Two Sample t-test
## 
## data:  MeanUptake[, i] by diagnose
## t = 1.6719, df = 18.286, p-value = 0.1116
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
##  -0.03026672  0.26745263
## sample estimates:
## mean in group HC mean in group PD 
##         2.147968         2.029375

## 
##  Welch Two Sample t-test
## 
## data:  MeanUptake[, i] by diagnose
## t = 2.6173, df = 21.182, p-value = 0.01603
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
##  0.02691038 0.23456181
## sample estimates:
## mean in group HC mean in group PD 
##         1.861839         1.731103

## 
##  Welch Two Sample t-test
## 
## data:  MeanUptake[, i] by diagnose
## t = 3.1029, df = 23.36, p-value = 0.004954
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
##  0.03958414 0.19752980
## sample estimates:
## mean in group HC mean in group PD 
##         2.051667         1.933110

## 
##  Welch Two Sample t-test
## 
## data:  MeanUptake[, i] by diagnose
## t = 3.2518, df = 23.352, p-value = 0.003468
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
##  0.03718395 0.16691831
## sample estimates:
## mean in group HC mean in group PD 
##         1.517768         1.415717

## 
##  Welch Two Sample t-test
## 
## data:  MeanUptake[, i] by diagnose
## t = -2.426, df = 17.77, p-value = 0.02615
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
##  -0.25449267 -0.01815847
## sample estimates:
## mean in group HC mean in group PD 
##         2.264076         2.400402

## 
##  Welch Two Sample t-test
## 
## data:  MeanUptake[, i] by diagnose
## t = -3.4362, df = 22.665, p-value = 0.002286
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
##  -0.26638508 -0.06607342
## sample estimates:
## mean in group HC mean in group PD 
##         1.909971         2.076201

## 
##  Welch Two Sample t-test
## 
## data:  MeanUptake[, i] by diagnose
## t = -2.8365, df = 24.157, p-value = 0.009087
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
##  -0.19255678 -0.03038962
## sample estimates:
## mean in group HC mean in group PD 
##         2.243670         2.355143
#association subcortical vs cortical left SN vs. cortical hypermetbolism precentral right
plot(gmFlexData$HC_gr_PD_m10m16m16,gmFlexData$HC_kl_PD_p5m28p61, col=gmFlexData$diagnose)
#cave: all data points 
model<-lm(gmFlexData$HC_kl_PD_p5m28p61 ~  gmFlexData$HC_gr_PD_m10m16m16)
abline(model)

ggsave("association_subcorticalvscortical.png")
## Saving 7 x 5 in image
summary(model)
## 
## Call:
## lm(formula = gmFlexData$HC_kl_PD_p5m28p61 ~ gmFlexData$HC_gr_PD_m10m16m16)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38537 -0.09142 -0.01741  0.06612  0.45476 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    2.41846    0.04547  53.182   <2e-16 ***
## gmFlexData$HC_gr_PD_m10m16m16 -0.29410    0.03147  -9.345   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1446 on 808 degrees of freedom
## Multiple R-squared:  0.09754,    Adjusted R-squared:  0.09642 
## F-statistic: 87.33 on 1 and 808 DF,  p-value: < 2.2e-16
#corr whole group mean values
cor.test(MeanUptake$HC_kl_PD_p5m28p61,MeanUptake$HC_gr_PD_m10m16m16)
## 
##  Pearson's product-moment correlation
## 
## data:  MeanUptake$HC_kl_PD_p5m28p61 and MeanUptake$HC_gr_PD_m10m16m16
## t = -3.3098, df = 25, p-value = 0.002836
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7703958 -0.2176128
## sample estimates:
##        cor 
## -0.5519803
#with ggplot
# Change point shapes and colors
ggplot(data=gmFlexData, aes(x=HC_gr_PD_m10m16m16, y=HC_kl_PD_p5m28p61)) +
  geom_point(aes(color=diagnose),size=3,alpha=0.5)+
    
  geom_point(data = MeanUptake, 
             aes(x = HC_gr_PD_m10m16m16, y=HC_kl_PD_p5m28p61, fill=diagnose),
shape=21, size=3, color="black")+
  geom_smooth(data = MeanUptake,
              method='lm', formula= y~x, color="black")+
   scale_color_manual(values=c("#1F78B4",  "#BFD200"))+
   scale_fill_manual(values=c("#1F78B4",  "#BFD200"))+
  theme_bw()+
     theme(legend.title = element_blank())+
  labs(title = "", x = "left SN",
           y = "Motor cortex")+
       theme(text = element_text(size = 30),
        axis.text.x = element_text(size=20),
        axis.text.y = element_text(size=20))+
   theme(legend.position = c(0.9, 0.8), legend.title=element_blank())+
    geom_text(x=1.6, y=2.47, label="r = -0.55, p = 0.003", color="black", size=10 )

ggsave("association_subcorticalvscortical_SN_M1.png")
## Saving 7 x 5 in image
#ROC curve
library(Epi)
## Warning: package 'Epi' was built under R version 4.3.3
MeanUptake$Fall<-MeanUptake$diagnose=="PD"
MeanUptake$Fall<-factor(MeanUptake$Fall)
str(MeanUptake)
## 'data.frame':    27 obs. of  15 variables:
##  $ subject            : Factor w/ 27 levels "RS_fPET001","RS_fPET002",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ HC_gr_PD_m10m16m16 : num  1.38 1.36 1.24 1.46 1.4 ...
##  $ HC_gr_PD_m14m98m4  : num  1.95 1.82 1.99 2.09 1.77 ...
##  $ HC_gr_PD_m16_m2_p12: num  1.78 1.7 1.28 1.74 1.63 ...
##  $ HC_gr_PD_m34_p48_p2: num  1.98 1.97 1.83 1.95 2.12 ...
##  $ HC_gr_PD_m52m58p24 : num  1.97 1.84 1.77 1.98 1.69 ...
##  $ HC_gr_PD_p14m96m16 : num  2.1 1.83 2 2.26 1.77 ...
##  $ HC_gr_PD_p22p2p10  : num  1.86 1.64 1.69 1.81 1.87 ...
##  $ HC_gr_PD_p46m62p36 : num  2.01 1.88 2 2.11 1.77 ...
##  $ HC_gr_PD_p8m16m16  : num  1.48 1.53 1.23 1.52 1.43 ...
##  $ HC_kl_PD_m33m9p4   : num  2.37 2.4 2.34 2.18 2.53 ...
##  $ HC_kl_PD_p5m28p61  : num  2.21 1.99 2.4 1.98 1.92 ...
##  $ HC_kl_PD_p35m2p3   : num  2.36 2.34 2.45 2.2 2.53 ...
##  $ diagnose           : chr  "PD" "PD" "PD" "HC" ...
##  $ Fall               : Factor w/ 2 levels "FALSE","TRUE": 2 2 2 1 2 1 2 2 2 1 ...
ResultsROC<-with(MeanUptake,
                 ROC(test=HC_kl_PD_p5m28p61, stat=Fall, plot="ROC",MI=F))

MeanUptake$HC_kl_PD_p5m28p61
##  [1] 2.209212 1.990884 2.401726 1.979455 1.918860 1.718265 2.010278 2.182069
##  [9] 2.128547 1.962481 1.963996 2.270269 1.839517 2.045838 2.003197 1.999362
## [17] 1.921927 1.897130 2.010274 1.994589 1.869792 1.723270 1.874417 2.054764
## [25] 1.940962 2.080473 1.904883
MeanUptake$Fall<-as.numeric(as.factor(MeanUptake$Fall))-1

#ROC curve with multiple variables
#all clusters
logit_model<-glm(Fall ~ HC_kl_PD_p5m28p61 + HC_gr_PD_p8m16m16+HC_gr_PD_m34_p48_p2+HC_kl_PD_m33m9p4+HC_gr_PD_m16_m2_p12+ HC_gr_PD_m10m16m16+HC_gr_PD_p14m96m16+HC_gr_PD_p46m62p36+HC_kl_PD_p35m2p3+HC_gr_PD_m14m98m4+HC_gr_PD_m52m58p24+HC_gr_PD_p22p2p10, data=MeanUptake)


MeanUptake$pred_probs<-predict(logit_model,type="response")
ResultsROC<-with(MeanUptake,
                 ROC(test=pred_probs, stat=Fall, plot="ROC",MI=F))
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

library(pROC)
## Warning: package 'pROC' was built under R version 4.3.3
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## 
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
#define object to plot
rocobj <- roc(MeanUptake$Fall,MeanUptake$pred_probs)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
auc <- round(auc(MeanUptake$Fall, MeanUptake$pred_probs),4) 
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
#create ROC plot
ggroc(rocobj, colour = '#00BFC4', size = 2) +
  ggtitle(paste0('ROC Curve ', '(AUC = ', auc, ')'))+
    theme_minimal()

#reduced clusters
logit_model_hypo<-glm(Fall ~ HC_gr_PD_p8m16m16+HC_gr_PD_m34_p48_p2+HC_gr_PD_m16_m2_p12+ HC_gr_PD_m10m16m16+HC_gr_PD_p14m96m16+HC_gr_PD_p46m62p36+HC_gr_PD_m14m98m4+HC_gr_PD_m52m58p24+HC_gr_PD_p22p2p10, data=MeanUptake)

MeanUptake$pred_probs_hypo<-predict(logit_model_hypo,type="response")

#increased clusters
logit_model_hyper<-glm(Fall ~ HC_kl_PD_p5m28p61+HC_kl_PD_p35m2p3+HC_kl_PD_m33m9p4, data=MeanUptake)

MeanUptake$pred_probs_hyper<-predict(logit_model_hyper,type="response")

#define object to plot
rocobj <- roc(MeanUptake$Fall,MeanUptake$pred_probs)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
rocobj2<-roc(MeanUptake$Fall,MeanUptake$pred_probs_hypo)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
rocobj3<-roc(MeanUptake$Fall,MeanUptake$pred_probs_hyper)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
auc <- round(auc(MeanUptake$Fall, MeanUptake$pred_probs),4) 
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
auc1 <- round(auc(MeanUptake$Fall, MeanUptake$pred_probs_hypo),4) 
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
auc2 <- round(auc(MeanUptake$Fall, MeanUptake$pred_probs_hyper),4) 
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
ci_roc1<-ci.auc(rocobj)
ci_roc2<-ci.auc(rocobj2)
ci_roc3<-ci.auc(rocobj3)

ci_band1<-ci.se(rocobj, specificities=seq(0,1,length.out=100))
ci_band2<-ci.se(rocobj2, specificities=seq(0,1,length.out=100))
ci_band3<-ci.se(rocobj3, specificities=seq(0,1,length.out=100))

ci_data<-data.frame(specificity=as.numeric(rownames(ci_band1)),
                    sensitivity=ci_band1[,2],
            lower=ci_band1[,1],
            upper=ci_band1[,3])

ci_data2<-data.frame(specificity=as.numeric(rownames(ci_band2)),
                    sensitivity=ci_band2[,2],
            lower=ci_band2[,1],
            upper=ci_band2[,3])

ci_data3<-data.frame(specificity=as.numeric(rownames(ci_band3)),
                    sensitivity=ci_band3[,2],
            lower=ci_band3[,1],
            upper=ci_band3[,3])

#create ROC plot
ggroc(list("All clusters" = rocobj,
           "Hypometabolic" = rocobj2,
           "Hypermetabolic" = rocobj3), size = 1) +
  ggtitle(paste0('ROC Curve ', '(AUC = ', auc, ')'))+
    theme_minimal()+
  
  geom_ribbon(data=ci_data,aes(x=specificity, ymin=lower, ymax=upper),
              inherit.aes= FALSE,fill="grey", alpha=0.2)+
  
    geom_ribbon(data=ci_data2,aes(x=specificity, ymin=lower, ymax=upper),
              inherit.aes= FALSE, fill=   "#1F78B4", alpha=0.2)+
  
   geom_ribbon(data=ci_data3,aes(x=specificity, ymin=lower, ymax=upper),
               inherit.aes= FALSE,fill= "#BFD200", alpha=0.2)+
  
   scale_colour_manual(values = c("All clusters" = "gray",
                                 "Hypometabolic" = "#1F78B4",
                                 "Hypermetabolic" = "#BFD200"),
                       guide = guide_legend(title = NULL))+

   theme(text = element_text(size = 30),
        axis.text.x = element_text(size=20),
        axis.text.y = element_text(size=20))+
  theme(legend.position = c(0.6, 0.4))

ggsave("roc_allClusters_ci.png")
## Saving 7 x 5 in image
setwd("~/Documents/RSfPETCog")

#load excel file
RSfPET<-read.xlsx(xlsxFile = "RS_fPET_fMRT_Data.xlsx",colNames = TRUE)

#define matrix for variation coefficient 
Vk<-matrix(ncol=27, nrow=12)

#loop over subjects to fill matrix with variataion coefficient per subject 

#calculate variation in subjects
Subjects<-c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029")


for(i in 1:length(Subjects)){
  
  Vk[1,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m10m16m16)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m10m16m16)*100
Vk[2,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m14m98m4)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m14m98m4)*100
Vk[3,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m16_m2_p12)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m16_m2_p12)*100
Vk[4,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m34_p48_p2)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m34_p48_p2)*100
  Vk[5,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m52m58p24)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m52m58p24)*100
  Vk[6,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p14m96m16)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p14m96m16)*100
  Vk[7,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p22p2p10)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p22p2p10
)*100
  Vk[8,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p46m62p36)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p46m62p36)*100
  Vk[9,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p8m16m16)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p8m16m16)*100
  Vk[10,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_kl_PD_m33m9p4)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_kl_PD_m33m9p4)*100
  
  Vk[11,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_kl_PD_p5m28p61)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_kl_PD_p5m28p61)*100
  Vk[12,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_kl_PD_p35m2p3)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_kl_PD_p35m2p3)*100
}

#transform into dataframe
Vk_data<-data.frame(
  Vk)
colnames(Vk_data)<-Subjects
tVk_data<-t(Vk_data)
colnames(tVk_data)<-colnames(gmFlexData[,1:12])
tVk_data
##            HC_gr_PD_m10m16m16 HC_gr_PD_m14m98m4 HC_gr_PD_m16_m2_p12
## RS_fPET001           9.149458          1.609397            5.520428
## RS_fPET002          11.174374          2.240504            8.320740
## RS_fPET003           8.022620          1.932647            9.629352
## RS_fPET004          10.962013          2.140445            7.782423
## RS_fPET005           8.504113          1.658767            7.300665
## RS_fPET006           8.744335          2.024273            8.351591
## RS_fPET007           8.111821          1.420442            4.606805
## RS_fPET008          11.383303          1.591634            8.671552
## RS_fPET010           9.868059          1.547594            7.268533
## RS_fPET011           5.203967          1.281799            5.507259
## RS_fPET012          10.621481          2.419876            7.434115
## RS_fPET014           9.893266          1.621636            7.851973
## RS_fPET015           7.897227          1.797752            6.324236
## RS_fPET016          11.171252          1.820831           10.433023
## RS_fPET017           6.924190          1.404396            6.623234
## RS_fPET018           8.918926          1.798407            7.407190
## RS_fPET019           8.683440          1.533015            6.715557
## RS_fPET020           6.211866          1.078035            4.450474
## RS_fPET021           8.110944          1.420442            4.606790
## RS_fPET022           8.671518          1.801752            5.561755
## RS_fPET023           6.272215          1.597641            5.153214
## RS_fPET024          14.587214          2.313922            9.904105
## RS_fPET025           8.964104          1.921321            7.927108
## RS_fPET026           9.633082          1.345150            7.113614
## RS_fPET027           7.512787          1.529448            8.411124
## RS_fPET028           7.130410          1.565766            8.438449
## RS_fPET029           7.721407          1.761468            7.009617
##            HC_gr_PD_m34_p48_p2 HC_gr_PD_m52m58p24 HC_gr_PD_p14m96m16
## RS_fPET001            2.116563           1.526527           2.986196
## RS_fPET002            2.339266           2.032020           2.874476
## RS_fPET003            1.428722           1.185980           2.478086
## RS_fPET004            2.413514           2.026804           2.747820
## RS_fPET005            2.327577           1.748098           2.923575
## RS_fPET006            2.092383           1.605132           2.224615
## RS_fPET007            1.266032           1.504451           1.456026
## RS_fPET008            2.676008           2.079021           2.604002
## RS_fPET010            2.194953           1.429761           2.496667
## RS_fPET011            2.107382           1.585552           1.429364
## RS_fPET012            6.745391           1.807977           2.899212
## RS_fPET014            1.734242           2.030424           3.190650
## RS_fPET015            1.974489           1.594484           2.193591
## RS_fPET016            2.329207           1.556435           3.368566
## RS_fPET017            1.796029           1.329104           1.613693
## RS_fPET018            1.611327           1.446840           2.177380
## RS_fPET019            1.571010           1.314334           1.852377
## RS_fPET020            1.162451           1.195506           1.225862
## RS_fPET021            1.265957           1.504512           1.456278
## RS_fPET022            1.672039           1.506442           1.868991
## RS_fPET023            1.973596           1.342139           2.538681
## RS_fPET024            2.003060           2.303311           3.850872
## RS_fPET025            1.761571           1.397084           2.770462
## RS_fPET026            1.919996           1.389658           2.613000
## RS_fPET027            2.469319           2.107585           2.203953
## RS_fPET028            3.513236           1.443132           2.615536
## RS_fPET029            1.713690           1.702189           2.463870
##            HC_gr_PD_p22p2p10 HC_gr_PD_p46m62p36 HC_gr_PD_p8m16m16
## RS_fPET001          5.171534           2.312951          7.675469
## RS_fPET002          6.612009           2.694771          9.149929
## RS_fPET003          5.026277           4.465347          9.592589
## RS_fPET004          6.786384           2.613275          9.908314
## RS_fPET005          5.199675           2.624397         12.276443
## RS_fPET006          6.292950           1.994300          9.857803
## RS_fPET007          4.210123           1.764034          5.636190
## RS_fPET008          5.625373           2.548775          8.576520
## RS_fPET010          5.609316           2.323284         11.026218
## RS_fPET011          3.923005           1.291514          6.123950
## RS_fPET012          6.284234           3.037558          8.272361
## RS_fPET014          6.907067           2.370053         10.131723
## RS_fPET015          5.269210           1.876998          9.049145
## RS_fPET016          6.026022           2.199640          9.985203
## RS_fPET017          4.178573           2.228051          6.765833
## RS_fPET018          4.344602           2.128189          9.290340
## RS_fPET019          4.015049           1.417233          8.160187
## RS_fPET020          4.014575           1.301726          7.246938
## RS_fPET021          4.210229           1.763910          5.636299
## RS_fPET022          4.971593           2.474279          8.798017
## RS_fPET023          4.479406           1.398938         10.423135
## RS_fPET024          8.101545           3.915617         15.915154
## RS_fPET025          4.355873           2.249852         11.326078
## RS_fPET026          5.901364           2.265688          7.772430
## RS_fPET027          5.116952           1.716365          9.864777
## RS_fPET028          4.337098           2.322549          9.663267
## RS_fPET029          4.831891           1.989895          8.642143
##            HC_kl_PD_m33m9p4 HC_kl_PD_p5m28p61 HC_kl_PD_p35m2p3
## RS_fPET001         6.605761         2.2677417         4.263425
## RS_fPET002         5.813709         2.0108309         4.140702
## RS_fPET003         6.637324         1.9997148         3.977227
## RS_fPET004         6.203789         1.5369858         4.011151
## RS_fPET005         5.917214         2.0229921         3.627611
## RS_fPET006         4.965329         1.3606797         2.920158
## RS_fPET007         4.702880         2.0652449         2.863745
## RS_fPET008         6.629785         2.2372681         3.851112
## RS_fPET010         4.722849         1.3858007         2.959791
## RS_fPET011         4.101901         0.8377055         2.757841
## RS_fPET012         8.266548         2.2360301         4.723090
## RS_fPET014         5.971227         2.5375296         3.942221
## RS_fPET015         5.650172         1.8496787         2.730661
## RS_fPET016         5.346292         1.5727147         4.449450
## RS_fPET017         4.666100         1.0056832         3.500288
## RS_fPET018         6.263513         1.7842069         4.216995
## RS_fPET019         5.991526         1.0462713         2.813564
## RS_fPET020         3.398688         0.9437292         2.754232
## RS_fPET021         4.703037         2.0654299         2.863897
## RS_fPET022         6.557544         1.3977717         4.405924
## RS_fPET023         5.684074         0.7907276         3.756343
## RS_fPET024         6.812369         1.6182024         7.303282
## RS_fPET025         5.584171         2.2873223         3.406044
## RS_fPET026         6.014291         1.1979635         3.045293
## RS_fPET027         6.803682         1.2740988         4.228673
## RS_fPET028         5.534481         1.3830695         2.675859
## RS_fPET029         4.566263         1.3834854         4.164847
#create group variable
#in groups HC vs PD 
ID_HCPD<-MeanUptake$diagnose

tVk_group<-as.data.frame(tVk_data)

tVk_group<-data.frame(lapply(tVk_group,as.numeric))
tVk_group$ID_HCPD<-ID_HCPD
tVk_group$ID_HCPD<-factor(tVk_group$ID_HCPD)
summary(tVk_group$ID_HCPD)
## HC PD 
## 13 14
variables<- colnames(tVk_group[,1:12])
grp<-levels(tVk_group$ID_HCPD)

#function
perm_test<- function(variable, group_var){
  group_var<-as.factor(group_var)
  grp_levels<-levels(group_var)
  
  obs_diff<-tapply(variable,group_var,mean, na.rm=TRUE)
  obs_diff_A<-obs_diff[grp_levels[2]]-obs_diff[grp_levels[1]]

  
  #perm
  perm_diffsA<-replicate(10000,{
    permuted_groups<-sample(group_var)
    perm_mean_diff_A<-tapply(variable,permuted_groups, mean, na.rm=TRUE)
    perm_mean_diff_A[grp_levels[2]]-perm_mean_diff_A[grp_levels[1]]
  })
  
  
  print(length(perm_diffsA))
  
  
  #p-values
  p_value_A<-mean(abs(perm_diffsA)>=abs(obs_diff_A))
  
  return(c(p_value_A))
  
}



#save results
results<-data.frame(Variable=character(), p_value=numeric(), stringsAsFactors = FALSE)

#loop over variables
for (var in variables){
  p_values<-perm_test(tVk_group[[var]], tVk_group$ID_HCPD)
  
  results<-rbind(results, data.frame(Variable=var, p_value_A=p_values[1]))
}
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
print(results)
##               Variable p_value_A
## 1   HC_gr_PD_m10m16m16    0.3454
## 2    HC_gr_PD_m14m98m4    0.9045
## 3  HC_gr_PD_m16_m2_p12    0.2495
## 4  HC_gr_PD_m34_p48_p2    0.8238
## 5   HC_gr_PD_m52m58p24    0.8775
## 6   HC_gr_PD_p14m96m16    0.0760
## 7    HC_gr_PD_p22p2p10    0.6586
## 8   HC_gr_PD_p46m62p36    0.1359
## 9    HC_gr_PD_p8m16m16    0.7014
## 10    HC_kl_PD_m33m9p4    0.7980
## 11   HC_kl_PD_p5m28p61    0.0235
## 12    HC_kl_PD_p35m2p3    0.9347
#gghalf plots vc 

library(gghalves)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# For multiple boxplots we need additional long format
long <- tVk_group %>%
  pivot_longer(
    cols = colnames(gmFlexData[, 1:12]), # Assumes gmFlexData has these columns
    names_to = "region",
    values_to = "value"
  )

# Calculate max value for each region to guide ordering
region_summary <- long %>%
  group_by(region) %>%
  summarise(max_y = max(value, na.rm = TRUE)) %>%
  arrange(max_y) # Arrange by ascending max value (regions with lower max_y first)
                 # Use arrange(desc(max_y)) for descending order if preferred

# Get the ordered region names
ordered_region_names <- region_summary$region

# Re-level the region factor in the long dataframe
# This new order will be respected by facet_wrap
long$region <- factor(long$region, levels = ordered_region_names)

# Create the plot
ggplot(long, aes(x = region, y = value, fill = ID_HCPD)) + 
  geom_half_violin(side = "r", trim = FALSE, alpha = 0.5, show.legend = FALSE) +
  geom_half_boxplot(side = "l") +
  facet_wrap(~region,scales="free", # You can add ncol or nrow here e.g. ncol=4 for 4 columns
             labeller = as_labeller(c(
               "HC_gr_PD_m10m16m16"  = "left SN", "HC_gr_PD_m14m98m4" = "left OC",
               "HC_gr_PD_m16_m2_p12" = "left AP","HC_gr_PD_m34_p48_p2" ="left MFC",
               "HC_gr_PD_m52m58p24"="left AG","HC_gr_PD_p14m96m16" ="right OC",
               "HC_gr_PD_p22p2p10"="right AP","HC_gr_PD_p46m62p36"="right AG",
               "HC_gr_PD_p8m16m16"="right SN" ,"HC_kl_PD_m33m9p4"="left PP",
               "HC_kl_PD_p5m28p61"="MC", "HC_kl_PD_p35m2p3"="right PP"
             )))+
   scale_fill_manual(labels=c("HC","PD"),values=c("#1F78B4",  "#BFD200"))+
  theme_minimal()+ 
  labs(title = "", x = "Region",
           y = "VC")+
  theme_minimal() + 
  labs(title = "", x = "Region", y = "VC") +
  theme(text = element_text(size = 25),
        title = element_text(size = 20),
        axis.text.x = element_blank(),
        axis.text.y = element_text(size=20),
        legend.title = element_blank())

# save plot
ggsave("Vk_datadriven_boxplots.png", width = 12, height = 10, dpi = 300)



#plot Vk vs clinical scores
tVk_data<-as.data.frame(tVk_data)
tVk_group$ID_HCPD<-factor(ID_HCPD)

#vs MoCA in whole group 
plot(tVk_group$HC_gr_PD_p14m96m16,subset(RSfPET,record_id!="RS_fPET013"&record_id!="RS_fPET009")$MoCa_Sum)

cor.test(tVk_group$HC_gr_PD_p14m96m16,subset(RSfPET,record_id!="RS_fPET013"&record_id!="RS_fPET009")$MoCa_Sum, method="spearman")
## Warning in cor.test.default(tVk_group$HC_gr_PD_p14m96m16, subset(RSfPET, :
## Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  tVk_group$HC_gr_PD_p14m96m16 and subset(RSfPET, record_id != "RS_fPET013" & record_id != "RS_fPET009")$MoCa_Sum
## S = 4615, p-value = 0.03428
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.408731
#   Spearman's rank correlation rho

#data:  tVk_data$HC_gr_PD_p14m96m16 and subset(RSfPET, record_id != "RS_fPET013" & record_id != "RS_fPET009")$MoCa_Sum
#S = 4615, p-value = 0.03428
#alternative hypothesis: true rho is not equal to 0
#sample estimates:
#      rho 
#-0.408731 


model<-lm(subset(RSfPET,record_id!="RS_fPET013"&record_id!="RS_fPET009")$MoCa_Sum ~  tVk_group$HC_gr_PD_p14m96m16)
abline(model)

summary(model)
## 
## Call:
## lm(formula = subset(RSfPET, record_id != "RS_fPET013" & record_id != 
##     "RS_fPET009")$MoCa_Sum ~ tVk_group$HC_gr_PD_p14m96m16)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.9036  -1.4685   0.6058   2.4538   4.2903 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    33.038      2.616  12.627 2.39e-12 ***
## tVk_group$HC_gr_PD_p14m96m16   -2.550      1.050  -2.429   0.0227 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.426 on 25 degrees of freedom
## Multiple R-squared:  0.1909, Adjusted R-squared:  0.1585 
## F-statistic: 5.899 on 1 and 25 DF,  p-value: 0.02268
#MC cluster 
cor.test(tVk_group$HC_kl_PD_p5m28p61,subset(RSfPET,record_id!="RS_fPET013"&record_id!="RS_fPET009")$MoCa_Sum, method="spearman")
## Warning in cor.test.default(tVk_group$HC_kl_PD_p5m28p61, subset(RSfPET, :
## Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  tVk_group$HC_kl_PD_p5m28p61 and subset(RSfPET, record_id != "RS_fPET013" & record_id != "RS_fPET009")$MoCa_Sum
## S = 4716.1, p-value = 0.02177
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4396019
#   Spearman's rank correlation rho

#data:  tVk_data$HC_kl_PD_p5m28p61 and subset(RSfPET, record_id != "RS_fPET013" & record_id != "RS_fPET009")$MoCa_Sum
#S = 4716.1, p-value = 0.02177
#alternative hypothesis: true rho is not equal to 0
#sample estimates:
#       rho 
#-0.4396019 

#with ggplot
ggplot(data=tVk_data,aes(x=HC_kl_PD_p5m28p61, y=subset(RSfPET,record_id!="RS_fPET013"&record_id!="RS_fPET009")$MoCa_Sum, color=ID_HCPD)) +
  geom_point()+
geom_smooth(method='lm', formula= y~x,show.legend = FALSE)+
  
    scale_color_manual(labels=c("HC","PD"),values=c("#1F78B4",  "#BFD200"))+
  theme_bw()+
   theme(legend.position = c(0.9, 0.8), legend.title=element_blank())+
  labs(title = "", x = "Motor Cortex VC",
           y = "MoCa sum score")+
   theme(text = element_text(size = 30),
          
        axis.text.x = element_text(size=20),
        axis.text.y = element_text(size=20))

#with Cog Z from demographic script
load("RSfPET.RData")

cor.test(subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$Cog_z,tVk_group$HC_kl_PD_p5m28p61, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  subset(RSfPET, record_id %in% c("RS_fPET001", "RS_fPET002", "RS_fPET003", "RS_fPET004", "RS_fPET005", "RS_fPET006", "RS_fPET007", "RS_fPET008", "RS_fPET010", "RS_fPET011", "RS_fPET012", "RS_fPET014", "RS_fPET015", "RS_fPET016", "RS_fPET017", "RS_fPET018", "RS_fPET019", "RS_fPET020", "RS_fPET021", "RS_fPET022", "RS_fPET023", "RS_fPET024", "RS_fPET025", "RS_fPET026", "RS_fPET027", "RS_fPET028", "RS_fPET029"))$Cog_z and tVk_group$HC_kl_PD_p5m28p61
## S = 4700, p-value = 0.02443
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4346764
#t = -2.7353, df = 25, p-value = 0.01129

#alternative hypothesis: true correlation is not equal to 0
#95 percent confidence interval:
# -0.7273049 -0.1222158
#sample estimates:
#       cor 
#-0.4799396 

#Figure
#with ggplot
ggplot(data=tVk_group,aes(x=HC_kl_PD_p5m28p61, y=subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$Cog_z,tVk_data$HC_kl_PD_p5m28p61, color=ID_HCPD, shape=RSfPET$MCI_NC_HC)) +
  geom_point(size=3)+
  
geom_smooth(data=tVk_data,aes(group=1),method='lm', formula= y~x,show.legend = FALSE, color="black", se=FALSE)+
  
    scale_color_manual(labels=c("HC","PD"),values=c("#1F78B4",  "#BFD200"))+
  theme_minimal()+
  # theme(legend.position = "none", #legend.title=element_blank())+
  labs(title = "", x = "Motor Cortex VC",
           y = "Cognition z-score",
       color="Group",
       shape="Group")+
   theme(text = element_text(size = 30),
        axis.text.x = element_text(size=20),
        axis.text.y = element_text(size=20))+
  scale_shape_manual(values=c("HC"=16, "MCI"=17, "PD-NC"=16,"PD-MCI"=17))+
  
geom_text(x=1.5,y=-1.0, label="r = -0.48, p = 0.011", color="black", size=7)
## Warning: The following aesthetics were dropped during statistical transformation: shape.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

 ggsave("Legend.png")
## Saving 7 x 5 in image
## Warning: The following aesthetics were dropped during statistical transformation: shape.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).